##*****************************************************************************
## FIGURE 3
##*****************************************************************************
## 
## R Version: R version 3.2.4 (2016-03-1)
##

##LOADING PACKAGES
library(ggplot2) ##Version: 2.1.0
library(MCMCpack) ##Version: 1.3-5 

##SET WORKING DIRECTORY
setwd('/workspace/fed_courts/scaling_paper/')

##
##LOADING DATASETS
##

##CERT DATABASE WITH ATTORNEY SCORES
cert.data <- read.csv('replication_materials/data/cert_petition_data_matched_to_dime.csv')


##WASH U SUPREME COURT DATA BASE, DOCKET CENTERED DATA
scdb <- read.csv("replication_materials/data/SCDB_2014_01_justiceCentered_Docket.csv")


##FILTERING CASES ON ROBERTS COURT  
scdb <- scdb[scdb$chief %in% c('Roberts') &
                 !(scdb$justice %in% c(104)) & ##EXCLUDING O'CONNOR VOTES (N=24, TOO FEW TO SCALE...)
                     scdb$term >= 2005,]
scdb$docket <- as.character(scdb$docket)

##
##MERGING SCDB AND CERT DATASETS
##
scdb.merged <- merge(x = scdb, y = cert.data, by = c("docket"),all=FALSE ,suffixes=c('','.y'))


##
##RECODING VOTE OUTCOMES -- (1 = VOTE FOR PETITIONER; 0 = VOTE FOR RESPONDENT)
##
scdb.merged$votePetitioner <- NA
scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 1] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 1] <- 0
 
scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 3] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 3] <- 0

scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 4] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 4] <- 0

scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 5] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 5] <- 0

scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 6] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 6] <- 0

scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 7] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 7] <- 0

scdb.merged$votePetitioner[scdb.merged$partyWinning == 0 & scdb.merged$vote == 2] <- 1
scdb.merged$votePetitioner[scdb.merged$partyWinning == 1 & scdb.merged$vote == 2] <- 0


##
##CONSTRUCTING VARIABLES WITH ATTORNEY IDEAL POINTS THAT CORRESPOND TO F
##
scdb.merged$winningcf.all <- scdb.merged$winningcf
scdb.merged$losingcf.all <- scdb.merged$losingcf

scdb.merged$winningcf <- NA
scdb.merged$winningcf[which(scdb.merged$votePetitioner == 1)] <- scdb.merged$cfscore_pet[which(scdb.merged$votePetitioner == 1)]
scdb.merged$winningcf[which(scdb.merged$votePetitioner == 0)] <- scdb.merged$cfscore_resp[which(scdb.merged$votePetitioner == 0)]

scdb.merged$losingcf <- NA
scdb.merged$losingcf[which(scdb.merged$votePetitioner == 0)] <- scdb.merged$cfscore_pet[which(scdb.merged$votePetitioner == 0)]
scdb.merged$losingcf[which(scdb.merged$votePetitioner == 1)] <- scdb.merged$cfscore_resp[which(scdb.merged$votePetitioner == 1)]

scdb.merged$cfpet <- NA
scdb.merged$cfpet[which(scdb.merged$votePetitioner == 1)] <- scdb.merged$cfscore_pet[which(scdb.merged$votePetitioner == 1)]
scdb.merged$cfresp <- NA
scdb.merged$cfresp[which(scdb.merged$votePetitioner == 0)] <- scdb.merged$cfscore_resp[which(scdb.merged$votePetitioner == 0)]


##
##CALCULATING AVERAGE SCORE OF PETITIONER ATTORNEYS BY JUSTICE
##
std <- function(x) sd(x)/sqrt(length(x))
rval <- NULL
for(ff in unique(scdb.merged$justice)){
    ##numv <- sum(!is.na(o1-o2))
    o11 <- scdb.merged$cfpet[scdb.merged$justice == ff]
    o12 <- scdb.merged$cfresp[scdb.merged$justice == ff]
    rval <- rbind(rval,c(justice.id=ff,
                         pet.mean=mean(o11,na.rm=T),
                         resp.mean=mean(o12,na.rm=T),
                         std.err.pet=std(na.omit(o11))
                         ##num.votes=numv
                         ))
}
rownames(rval) <- scdb.merged$justiceName[match(rval[,'justice.id'],scdb.merged$justice)]


###############################################################################
##FIT MCMC IRT MODEL
###############################################################################

##FILTERING ON CASES WHERE IDEAL POINTS ARE AVAILABLE FOR BOTH ATTORNEYS 
scd <- scdb[scdb$docket %in% scdb.merged[!is.na(scdb.merged$winningcf) & !is.na(scdb.merged$losingcf) ,'docket'],]


##RESHAPING SCDB DATA INTO VOTE MATRIX 
ujn <- unique(scd$justiceName)
ucase <- unique(scd$docket)
mqrc <- matrix(NA,length(ucase),length(ujn))
colnames(mqrc) <- ujn
rownames(mqrc) <- ucase
for(ff in unique(scd$justice)){
    dd <- scd[scd$justice == ff,]
    m1 <- match(rownames(mqrc),dd$docket)
    m2 <- match(unique(dd$justiceName),colnames(mqrc))
    mqrc[,m2] <- dd$vote[m1]
}

##RECODING VOTE OUTCOMES
mqrc[mqrc==3] <- 1
mqrc[mqrc==4] <- 1
mqrc[mqrc==5] <- 1
mqrc[mqrc==6] <- 2
mqrc[mqrc==7] <- 2
mqrc[mqrc==2] <- 0
mqrc[!(mqrc %in% c(0,1))] <- NA

##MCMC IRT SCALING MODEL--20,000 BURN-IN, 100,000 ITERATIONS
posterior <- MCMCirt1d(t(mqrc),
                        theta.constraints=list(JPStevens="-",
                            AScalia="+"),
                       seed=12345,
                       burnin=20000,
                       mcmc=100000,
                       store.item=TRUE,
                       thin=10, verbose=500)

##EXTRACTING ESTIMATED IDEAL POINTS FOR JUSTICES
mq <- t(rbind((colMeans(posterior)),(apply(posterior,2,quantile))))
votes <- mq[grep('beta',rownames(mq)),]
rownames(votes) <- gsub('beta\\.','',rownames(votes))
mq <- mq[grep('theta',rownames(mq)),]
rownames(mq) <- gsub('theta\\.','',rownames(mq))
mq <- mq[match(rownames(rval),rownames(mq)),] 
mq[,4] <- as.numeric(mq[,4])
mq <- as.data.frame(mq)

##SCATTERPLOT
pet.atty.mean = rval[,'pet.mean'] 
q <- qplot(x=pet.atty.mean,y=as.numeric(mq[,4]),label=rownames(mq))
q <- q + geom_text(aes(x=pet.atty.mean,y=as.numeric(mq[,4])+.06))
q <- q + geom_point()
q <- q + xlab('Avg. Ideal Point of Petitioner Attorneys Sided With')
q <- q + ylab('Martin-Quinn IRT')
q <- q + xlim(min(pet.atty.mean)-.05,max(pet.atty.mean)+.05) + ylim(min(mq[,4])-.065,max(mq[,4])+.065)
q <- q +  theme_bw()
q <- q + theme(axis.title= element_text(size=14))

##WRITING PDF TO DISK
pdf(file='replication_materials/figures/figure_3_avg_atty_ip_v_mq_irt.pdf',width=9,height=9)
print(q)
dev.off()

##WRITING EPS TO DISK
setEPS()
postscript(file='replication_materials/figures/figure_3_avg_atty_ip_v_mq_irt.eps',width=9,height=9)
print(q)
dev.off()



